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DETERMINING VISIBILITY INTERVALS BETWEEN AN EARTH 
TRACKING STATION AND A PLANETARY SATELLITE 


SUMMARY 


The problem of finding visibility intervals between an Earth 
tracking station and a planetary satellite is formulated based on the 
assumption of an infinite speed of light. A modification to account 
for the finite speed of light is then added to the formulation. An 
input-output description and listing of a Fortran computer program 
which calculates visibility intervals is provided in Appendix 1. The 
program was used to obtain visibility intervals between the Gold- 
stone deep space tracking station and a satellite in a typical elliptic 
orbit around Mars. 




DETERMINING VISIBILITY INTERVALS BETWEEN AN EARTH 
TRACKING STATION AND A PLANETARY SATELLITE 


INTRODUCTION 

The purpose of this paper is to provide an efficient technique for determining 
intervals in which a planetary satellite is visible to a given earth tracking station. 
This problem is similar to the one treated by Stern (1) who determines communi- 
cation times between an earth satellite and a planetary satellite. There is also 
a similarity to the problem of finding shadow times of a satellite-Yeremenko (2) 
Karytevn (3), and others. But the problem treated in the present paper is not a 
special case of the ones mentioned above and hence it requires its own separate 
formulation. 

Clearly the visibility between an earth tracking station and a planetary satel- 
lite can be obstructed by both the earth and the planet. The approach taken will 
be to first derive an earth occultation function, S £ , of time from epoch which is 
positive for times at which the earth fails to occult visibility and negative for 
times at which the earth does occult visibility. A similar function S p for plane- 
tary occultation is also derived. A visibility function V of time from epoch can 
then be defined which is equal to one when S £ and S p are positive, and is equal 
to minus one otherwise. The discontinuities of V can be easily and accurate^ 
determined by numerical methods. The times associated with these discontinuities 
are either visibility rise or visibility set times depending on the nature of the 
discontinuity. By this technique, visibility intervals or "windows" can be constructed. 

The earth and planetary occultation functions will be derived assuming an 
infinite speed of light. This permits the application of direct line of sight geometry 
in the derivations. Later a slight modification of the occultation functions will 
substantially correct for the error introduced by this assumption. 

In Appendix I is an input and output description and listing of a Fortran 
routine called "OCCULT" which determines visibility windows between an earth 
tracking station and a planetary satellite. The program is based on the analysis 
outlined above. This program will be utilized in a later section to determine 
visibility windows between an earth tracking station and an artificial satellite in 
a typical two body orbit around Mars. 
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DERIVATION OF EARTH OCCULT ATION FUNCTION 

The object of this section is to define a function S £ of time which assumes 
negative values for times when the Earth occults visibility between an Earth 
tracking station and a satellite in orbit around a planet, and which assumes positive 
values otherwise. The derivation will be carried out under the following 
assumptions: 

(1) The speed of light is infinite. 

(2) The satellite is in a two body elliptic orbit about the planet. 

(3) Earth nutation can be ignored. 

(4) The Earth's angular velocity is constant. 

(5) The Earth and the planet are perfect spheres. 

(6) The interval in which occultation information is desired is small enough 
so that the position vector of the planet relative to the Earth can be con- 
sidered independent of time. 

(7) All vector quantities are given with reference to a mean equinox of date, 
mean equator of date coordinate set where the epoch date is a given time 
of perifocal passage of the planetary satellite. 

(8) Along with latitude and longitude, another parameter associated with a 
given tracking station is its elevation angle 8 • Visibility between a 
tracking station and a satellite will be considered occulted when the 
angle between the tangent plane at the tracking station and the line of 
sight to the satellite is less than b . The angle 8 is typically around 5 
degrees. 

Figure 1 displays the geometric situation at the exact instant of Earth 
occultation. The position vector from Earth center to tracking station is repre- 
sented by R £ . The vector R is the position vector of the planet relative to the 
Earth center. R is the position vector of the satellite relative to the planet 
center and R is the position vector of the satellite relative to the Earth center. 

The dotted line shown in Figure 1 is the line of sight from the satellite to the 
tracking station. The angle between this line and the tangent plane at the tracking 
station is the elevation angle 8 . The symbol y represents the angle between the 
position vectors of the tracking station and satellite relative to Earth center and 
at all times is given by 
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Figure 1 -Geometry for Earth Occupation Function 



cos H* - Rg • R s /| R E | | R, 


(1) 


where 


R g - R + Rp 


(2) 


At the instant of Earth occultation V satisfies the relation 


cos ( S + ¥) = j R | cos 5 / | R, 


(3) 


With the aid of some trigonometric identities, equation 3 can be written as 


cos 8 cos V - sin 8 yi - cos J, P - | R | cos 8 ' | R_ | = 0 


(4) 


By substituting the right side of 1 for cos y in 4, the Earth occultation function 
is obtained 


cos 8 ( R £ ■ R s ) - sin 8 R £ l 2 I R s 1 2 - ( R E ' R s ^ 2 - I R £ I 2 cos 8 ^ ^ 

Ir e I I R s I 

S £ is clearly zero at Earth occultation. When the satellite is visible to the track- 
ing station, S £ is positive as can be seen by examining 5 when V = 0, a condition 
which clearly implies visibility. During Earth occultation S £ is negative as can 
be seen by examining Equation (5) when V = n /2. 

Substantial computing time can be saved by noticing that when is between 
fl/2 and 3n/2, visibility is occulted by the Earth. Hence when r e - R g < 0, 

can be assigned some arbitrary negative value without altering the desired 
characteristics of S £ . The adjusted definition of the Earth occultation function 

S E is 


S 


E 


cos S (P £ • R s ) - s in 8 


^| R £ | 2 | R s | 2 - (R e R s ) 2 - | R £ | 2 cos 5 

R £ I |Rg 


( 6 ) 


when R £ R s >0 and S £ = - 1 otherwise. 
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Equations (1) through (6) recursively define S £ as a function of R E and R p . 

In order to define S £ as a function only of time from an epoch date it is necessary 
to define R £ and R p as functions of time from an epoch date. 

Let H be the epoch date, where M is a time of periapsis passage of the plane- 
tary satellite. Lot T be days from epoch. Let LAT and LONG be the latitude and 
longitude of the tracking station in radians. Also let H be the right ascension of 
Greenwich meridian from mean equinox of epoch time 3. Then R £ is given by 


R £ ( 1 ) t R n cos (LAT) cos (LONG + H + 2llT) 


R £ ( 2) = R 0 cos (LAT) sin (LONG + H + 2HT) 

R £ ( 3) - R fl sin (LAT) 
where R 0 is the radius of the earth. 

To obtain R p as a function of time from epoch, the orbital elements of the 
satellite must be available. Let a> - argument of perigee of satellite orbit, 0 = 
longitude of ascending node of satellite orbit, i = inclination of satellite orbit, 
a - semi -major axis of satellite orbit, c = eccentricity of satellite orbit. These 
orbital elements are assumed to be given relative to a mean equinox of epoch 
mean equator of epoch coordinate set. The unit vectors p and 0 are defined as 


P(1 ) - cos & cos Q — sin co sin Cl cos i 
P( 2 ) = cos uj sin fi + sin ^ cos fleos i 
P ( 3 ) - sin c*j sin i 


Q(l)-_sinc*j cos Q - cos w s i n fi cos i 

Q( 2 ) = - s i n 6 l> s in Q + cos a> cos Q cos i 
Q( 3 ) - - cos oj sin i 

In terms of eccentric anomaly E of the satellite, R p can be represented as 

R p - a( cos E + e ) P + a( 1 - e 2 )* 2 s in E 0 (9) 
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The functional relationship between eccentric anomaly E and time from epoch T 
is given by 


nT - E - o sin E 


( 10 ) 


where n is the mean motion of the satellite orbit. Equation (10) transendentally 
defines E in terms ofT . Hence to obtain the corresponding value of E for a 
given value of T requires a numerical process. 

Equations (8), (9), and (10) recursively define R p as a function of days from 
epoch T. Thus Equations (1) through (10) recursively define S p as a function of 
T where T represents days from an epoch time T, a time of satellite pcriapsis 
passage. The function S E (T) is zero at the instant of initiation or termination of 
Earth occultation. S E (T) is negative for times in which Earth occultation occurs 
and is positive for all other times. 


DERIVATION OF PLANET OCCULTATION FUNCTION 

In this section a function S p of time from epoch will be derived which is 
negative when the line of sight between a planetary satellite and an Earth tracking 
station is occulted by the planet and is positive otherwise. The same assumptions 
listed in the derivation of the Earth occultation function will be in force here. 

Figure 2 displays the geometric situation at the instant of initiation or termi- 
nation of planet occultation. The dotted line is the direct line of sight between the 
satellite and the tracking station. R E is the position vector of the tracking station 
relative to Earth center. R represents the position vector of the planet relative 
to Earth center. R t is the position vector of the tracking station relative to the 
planet center. R p is the position vector of the satellite relative to the planet. 

The symbol V represents the angle between R ( and R p . is the acute angle 
between R t and the line of sight. <t> is the angle between R p and the line of sight. 
At all times ^ is given by 


cos Y - R t • R p ' lR t | |R p | 


(ID 


where 


R t -R e -R 


( 12 ) 


At an instant of initiation or termination of planetaij occultation, V satisfies the 
relation 
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q» = n - Q - 4> 


(13) 


where 9 <n /2, ^ < n /2 and 


cos 0 = x p / I R t I , cos 0 T, x p / | R p I 


(14) 


where X p is the radius of the planet. These facts together with certain trigonometric 
relations yield 


cos ¥ - 


V(|RJ 2 -X p 2 )( R | - X 2 ) - X 2 


l R J i R p l 


(15) 


The difference between the right sides of Equations (11) and (15) is the planetary 
occultation function 


S = 

p 


R • R + X 2 - 

t p p 


V(|R t i 2 -X 2 )(|R p | 2 -X 2 ) 


(16) 


IRJ R. 


At an instant of initiation or termination of planetary occultation S p is zero. When 
the planet is not occulting viaibility, s is positive as can be seen by examining 
Equation (16) when ¥ = 0. When the planet occults visibility S p is negative as 
can be seen by examining Equation (16) when =n . 

Computational time can be saved by noticing that when -11/2 < y < n /2, 
planetary occultation cannot occur. Thus the planetary occultation function can 
be redefined as 


S = 

p 


_ R t • Kp + \ 2 - I R t l 2 - X 2 ) ( |R p | 2 - X 2 ) 


(17) 


IRJ IRpl 


when R • R <0 and S =1 otherwise. 

t p p 


4 


8 


Equation (17) along with Equation (12) recursively define S p as a function of 
R and R_. . But R_ and R can be expressed as functions of time from epoch as 

PE E p 

was shown in the previous section. Hence the derivation is complete. 


ON THE CALCULATION OF VISIBILITY INTERVALS 

The possession of Earth occultation and planetary occultation functions makes 
the calculation of visibility intervals an easy matter. A function v can be defined 
in the following way 

V ( T) = 1 if S_(T) >0 and S „(T) > 0 

£ p 

V(T) = - 1 otherwise 


The intervals in which V equals one correspond to visibility intervals. The 
calculation of these intervals is just a question of discovering the location and 
type of the discontinuities of V . This can be done to any desired accuracy by the 
application of simple numerical procedures. 

The independent variable in the visibility function V is time in days from an 
epoch where the epoch time is a time of periapsis passage of the planetary satel- 
lite. The visibility function is not defined for negative values of T . A uniformly 
effective procedure for the calculation of visibility intervals within a given time 
interval is now evident. The time interval i; which occultation information is re- 
quested must start with a time of periapsis passage 3 . The interval ends at some 
later time 3 + T Q . To obtain visibility intervals during this period it is necessary 
to obtain the discontinuities of the visibility function V in the interval [ 0, T Q ] . 

A straightforward numerical procedure incorporated in a computer program can 
determine these discontinuities and also determine their type. If values im- 
mediately to the left of a discontinuity are negative, then the time associated with 
the discontinuity is the beginning of a visibility interval. Otherwise a time 
associated with a discontinuity represents the end of a visibility interval. 

Some of the assumptions made in the derivation of the Earth and planetary 
occultation functions place restrictions on the accuracy and effectiveness of the 
procedure outlined above. Notice that no provisions are made for updating the 
position vector of the planet relative to the Earth during the interval in which 
visibility information is requested. Thus the longer this interval is, the more 
serious the error becomes. Obviously this presents no serious difficulty since 
any interval in which visibility information is required can be divided into inter- 
vals small enough so that this error can be ignored. The position vector of the 
planet relative to the Earth would then be updated at the beginning of each of these 
intervals. 


Another source of inaccuracy is the assumption of an infinite speed of light. 

If the distance between the Earth and the body about which the satellite is orbiting 
is great then the error caused by this assumption can be significant. But a slight 
alteration in the definitions of the Earth and planetary occultation functions can 
substantially correct for this error. This is the subject of the next section. 


A MODIFICATION TO ACCOUNT FOR THE FINITE SPEED OF LIGHT 

The assumption of an infinite speed of light in the derivation of the occultation 
functions permitted the application of line of sight geometry. But it introduces 
an error which can be significant when the Earth and the planet in question are 
separated by a great distance. Stern (1) in determining communication times 
between an Earth satellite and a planetary satellite effectively deals with this 
problem. 

The treatment provided here is simpler than that given by Stern. It rests 
on the assumption that the distance between the satellite and the tracking station 
can be well approximated by the distance between the Earth and the planet. Thus 
the time of travel in days of the signal from the satellite to the tracking station 
is AT= |r|/C where again R is the position vector of the planet relative to the 
Earth and C is the speed of light in kilometers per day. To determine if occulta- 
tion occurs at a time T , line of sight geometry can again be used provided the 
line of sight is understood to be between the tracking station at time T and the 
satellite at time T - AT. Hence in order to correct the occultation functions so 
as to account for the finite speed of light, all that is necessary is to compute 
the R p vector of Equations (2) and (17) at time T - AT instead of time T where 
T is the input variable of the occultation functions. 
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AN EXAMPLE 

In Appendix 1 is an input and output description and listing of a Fortran 
computer program called "OCCULT.” This program calculates visibility intervals 
between an Earth tracking station and a planetary satellite. The calculation is 
performed by obtaining the visibility function described in previous sections and 
locating its discontinuities. These discontinuities are then properly interpreted 
as either initiation or termination points of visibility. From this information, 
visibility intervals are then constructed and provided as output. 


This program was used to obtain visibility intervals between the Goldstone 
deep space tracking station and a satellite in a typical Mars orbit. In order to 
make the geometric situation as clear as possible, Mars was assumed to lie on 
the principal axis of the mean equinox mean equator of epoch coordinate set. 

This, of course, makes the example an artificial one, but it permits the reader to 
grasp the relationship between the Earth-Mars line and the orbit of the satellite 
since the orbital elements ai e also given with respect to the mean equinox mean 
equator of epoch coordinate set. Mars was assumed to be two hundred million 
kilometers from Earth and the satellite orbit was chosen as to insure that planetary 
occupation would occur. The parameters involved in the calculation are given below. 


Planet — Mars 

Radius = 3393.4 KM 

Gravitational constant = 42977.8 KM 3 /sec 2 
distance from Earth = 200,000,000 KM 

Tracking Station - Goldstone 

longitude = 243.11 degrees 
latitude = 35.24 degrees 
occupation angle = 5 degrees 


Satellite Orbit - Mars Elliptic 


Semimajor axis = 5 Mars radii 

eccentricity = .5 

inclination = 5 degrees 

longitude of ascending node = 5 degrees 

argument of perigee = 45 degrees 


Interval in which occupation information is requested = 3 days. 


| 




The results obtained by the OCCULT program are displayed graphically in 
Figure 3. The solid rectangles on the bottom scale illustrate the time intervals 
in which visibility is permitted between the Goldstone tracking station and the 
Mars satellite. The solid rectangles shown on the middle scale indicate the 
intervals in which the Earth does not interfere with visibility between tracking 
station and satellite. The solid rectangles on the top scale 'ndicate the intervals 
in which Mars permits visibility between tracking station and satellite. As one 
would expect, Earth occultation accounts for the major portion of lost visibility 
time between tracking station and satellite. In this particular case, during the 
three day period under investigation Earth occultation accounted for over 93% 
of lost visibility time. It also should be mentioned that the OCCULT program 
which performed these calculations contains a modification to account for 
the finite speed of light. 


CONCLUSION 

An effective procedure for determining visibility intervals between an Earth 
tracking station and a planetary satellite has been given in this report. The pro- 
cedure has been incorporated into a computer program, and input and output 
description and listing of which is given in Appendix 1. The assumption of an 
infinite speed of light is a convenience but it can in some situations lead to a 
significant error. A modification to account for the finite speed of light is 
easily performed. 

The program listed in /">Dendix 1 was used to calculate visibility intervals 
between the Goldstone deep - >ace tracking station and a satellite in a typical two 
body orbit around Mars. The results are displayed graphically in Figure 3. 

In this particular case, Earth occultation accounted for over 93% of lost visibility 
time during the three day period under investigation. 
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APPENDIX 1 

Input-CXitput Description and Listing of "OCCULT" 


PROGRAM DESCRIPTION OF "OCCULT" 

OCCULT is a double precision Fortran program written for the IBM 360 
model 91 or 95 computer. Its purpose is to compute visibility windows between 
an earth tracking station and a probe executing a closed two body orbit around 
another body. The program assumed a finite speed of light. An interpretation 
of the arguments in the calling sequence of OCCULT follows. 

CALL OCCULT (EPS, A, XI, SIGMA, OMEGA, U, XP, RP, 

XIT, TI, TEND, LAT, L0NG, DELTA, YT, T, N, OPEN) 


INPUT 


EPS - Eccentricity of satellite orbit about planet 

A - Semimajor axis of satellite orbit about planet (in KM) 

XI - Inclination of satellite orbit about planet (in radians) 

SIGMA - Longitude of the ascending node of satellite orbit about planet (in 
radians) 

OMEGA - Argument of Perigee of satellite orbit about planet (in radians) 

U - Gravitational constant of planet in (KM 3 /sec 2 ) 

XP - Radius of planet (in KM) 

RP - A three dimensional array giving the rectangular coordinates of 

the planet relative to a mean equinox mean equator of epoch coordinate 
set (in KM) 

XIT - Integer number of days from January 1, 1950 to a time of perifocal 
passage 

TI - Fractional number of days from January 1, 1950 to a time of perifocal 
passage 

TEND - Number of days after perifocal passage in which occultation ifnormation 
is requested 
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LAT 

LONG 


DELTA 


- Latitude of tracking station (in radians) 

- Longitude of tracking station (in radians) 

- Elevation angle of tracking station (in radians) 

- Time interval, in days, smaller than any visibility window of interest 

- Accuracy desired, in days, in visibility window times 


OUTPUT 

j = -1 if there are no visibility windows in interval 
N = 0 if there is no occultation in interval 

L= Number of visibility windows in interval otherwise 

OPEN -An N by two array: 

OPEN (I, 1) is rise time in days from perifocal passage of I th 
visibility window 

OPEN (I, 2) is set time in days from perifocal passage of I th 
visibility window 

IfN = 0orN = -l, the OPEN array is undefined. 

The OPEN array lists the visibility windows in chronological order. 
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X R P ( T sRP ( I ) 
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ARPS=DSORT ( ARPS ) 


i in. m 


ROHTsRS ( 1 )*RPS ( 1 ) +RS( 2 )-KPS( 2 )+RS( 3) *RPS( 3) 


n . nnn i r.i i t 1 1 



GO TO 10 

6 7_ = r)S'-.n<T ( ( AKS*ARS-XP*XP )» ( AKPS*Ak 


$? = (RQOT+Z-(XP*XP ) ) /( AR SHARPS) 

0 wr 
F Ain 

R f A L FUNCTION YRUOT*P 


MXllIVlilfl 


tart a stop of intfkva 


T = aCCURACY DESIRED, F UM C = FUNCTION JAwF TO DEFINE 




iwiHKiHBaanH 


ilKIHWia 



OUTPUT VARIABLES I OUT* NUMBER OF ROOTS F0UMQ ( 'aX=100) 

Z E R 0 = A R R AY OF ROUTS ( m AX = li»0 ) * aY = AR R A Y uF RIX. 


INDICATORS WHERE +1.0DO=SET AND -1.00U=RISE ROOT 

I iPLICIT REAL * 8 (a-H,0-Z) 


n I m E imS I On ZERO (100), WAY (100) 
DUMr Y=99.0D0 


X N = M 

M = WAX NUNBER OF ITERA 


) - P ( 1 ) 
E.XP) GO TO 


•.= ( OLOG(XP) -DLOG(XM) -OLDG(T) ) /0L0G( 2.000) 
IF 


VALUE = P(l) 

I ii IT = Q 


k I GHT=+ 1 .000 

L E FT = R I GOT 


TEST* FUNC (VALUE ) 

— OCOS FUNC WAS USED FOR CHECK-OUT Pi'RPQS 


IF ( TEST. GE. 0.000) GO TO 9 

L FFT=- 1 . 0D0 


R TGHT = L EFT 

OFLTA=(P (2)-P( 1) ) /XM 


\/ALOE = p ( i ) -DELTA 

'■>0 10 ISTEP = 1 . !'■< 


VALOE=VALUE + delta 

KFE P = L E FT 


T C ST = FUNC (VALUE) 
R T GUT = + 


IF ( TEST .LT .0.000) RIGHT* -1.000 

IF (LEFT. Me. RIGHT) GO TO 11 

G n TO 11 . . . CURVE HAS CROSSED AXIS 

GO TO 10 










11 XMARK=VALUE 
STOP* VALUEL 
CENTER* ( VALU 


icsiiajg 

RVmi 


IF(TEST.LT.O.ODO) SIGN=-1.000 

i 


ESfld 


; ; 


CENTER* ( CENTER - 
on TO 1 4 

XMARK) /2,0D0 + XMARK 

i 

12 

XMARK=CENTER 


■ 




* 

14 

C 0 N T I N U E 


/ 

100 



ZER0{ I0UT )*CENTER 


^VIUlUlgH» 

affinl 


LEFT*+1.0D0 


RIGHT=LEFT 






GO TO 778 
R 


|4aanMHianpa»iir^ini«g 

luMsaWI 

IB 

uiam 


FUNCTION E A ( E ) 

* 


REAL * 8 LONG »L AT 


,XT 


EA*E-(XN*YTIME)-(EPS*DSIN(E) ) 
R 


END 
SUB 


IMPLICIT REAL * 8 <A-H,0-Z) 


D=T'/J 

T=TF*864Q0 


0MEGA*0 ,72921 15070-4 - 0.38D-16 * D 


X 0.50641D- 14), 6, 283185307 17958600) 
RETURN 





F< MvCT I ON V 1 Z ( T I - E 


ILLICIT REAL * 8 (A-H,0-Z) 

REAL * 8 LONG *L AT 


C • V , - on GHAI ,RX,L ON G,LAT, XW , FPS, A, P, 0, DELTA, P.P, XP 
1 » X T 


0 I N FNS I (Jim R ( 3 ) , RPS (3) , R S ( 3 ) , P ( 3 ) , 0( 3 ) , R P ( 3 ) 

YTIME=TImE-XT 

IF(S1 (TIME) .LT. O.ODO) GU TO 10 

IF(SZ(TI.mE) .LT. Q.QDO) £U_T D .1,0 

'/ 17 = 1.000 

GO TO 15 


10 R I Z = - 1 • ODD 

18 R ETURN . 

I- h i " 

REA L FUNCTION ROOT»8 (N,P« T « TOUT, FUMCt ZERO. PAY ) 

C- FUNCTION S/R TO COMPUTE THE ZERO CROSSINGS OF A 

C CONTINUOUS CURVE ON CART COORDINATES 

C CAM BE USED AS A FUNCTION OR A S/R 

C INPUT VARIABLES N=NUMRER OF STEPS P=TR*0 PLACE ARRAY 

C WITH START £ STOP OF INTERVAL TO BE EXAMINED 

C T=AGCUR AC Y DESIRED, FUNC= FUNCTION MAoE TO DEFINE 

C- CURVE (PASSED W/ EXTERNAL STATEMENT IN CALLING PROG) 

C OUTPUT VARIABLES IOUT=WMRER OF ROOTS FOUND ._( MAX= 10Q ) , 

C ZERQ=ARRAY OF ROOT S ( HAX = 100 ) , WAY= ARRAY OF ROOT RISE/SFT 

C INDICATORS WHERE +1.0P0=SFT AND -1.0D0=RISE ROOT 


IMPLICIT REAL * 8 ( A-H , 0-Z ) 

DIMENSION-ZERO (100), WAY _ LIQDJ . _P_( 2J 

DUKMY=99.0D0 

X M = M 

C M = MAX NUMBER OF ITERATIONS AMD IS A FUNC OF P,N,T 

XP= P< 2) - P(1I 

IF(T.GE.XP) GO TO 777 

M= ( DL OG ( X P ) -PLUG ( XN ) -PLUG ( T) ) /DLOG( 2. ODD ) 

IF(M.LE.O) GO TO 777 

VA LUE > 

inuT=o 

R I G HT = + L.D DQ _»_ ■ 1 

LFFT=RIGHT 

TE ST=FUN C..( V ALU E J 

C --DC OS FUNC WAS USED FOR CHECK-OUT PURPOSES 

IFtTEST .»EE.iQLt O Qfl I SO T O -9 - 

LFFT=-1.0D0 

RIGHTfLEFT 








HI 


nElTA«(P(Z)-P(l))/XN 
VALUE®P ( 1 ) -DELTA 


DO 10 ISTEP* 1 ,N 
VALUE® VaL UE + DELTA 


KE 6 P=LE FT 
TFST=FUNC (VALUE) 


IF(lEFT.NE. RIGHT) GO TO 11 
Gn TO 11 ... . CURVE HAS 


GO TO 10 

XMARK® VALUE - DELTA 


DO 14 I T E R = 1 » H 
TEST- FUNG (CENTER) 


ODO) SIGN® 


IF(SIGM.EO.LEFT) GO TO 12 




DO + X 1 ' -ARK 


VALUE=C ENTER 
CENTER® ( CENTER - X 


Gn TO 14 
XMARKsCENTER 


CENTER ® ( VALUE - CENTER) / 2.000 4 CENTER 
CUNT I N U E .. . . . . 


I OUT = I OUT + 1 
ZERO! 1 QUT)=CENTER 


N A Y ( I OUT) =KEFP 
L 6 FT®+ 1 *QDG 


ROOT = DUNN Y 


ZERO ( 1 ) = < P ( 2 ) -P ( 1 ) )/ 2 .QD 0 

-1 


R DOT® DUNN Y 
RETURN 


m 
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liwlisf 


mm 


pG.j 

mm 


H 


■H 






jama 


HR 
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